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Computation  of  Supersonic  Space  Encircling  Plow  of  Blunt-Nosed  Body 
Zhu  Youlan,  Wang  Ruquan,  Zhong  Xichang 

Computer  Technology  Research  Institute 
The  Chinese  Academy  of  Sciences 

I.  Introduction 

Since  the  1950's,  "or  numerical  solution  of  the  problems  of  supersonic 

encircling 

inviscid  ^  flow  of  blunt-nosed  body,  a  number  of  different  methods 

has  been  developed.  Of  them  one  categorv  is  stationary  method  and  the 
other  is  nonstationary  method.  In  the  category  of  stationary  method,  there 
are  method  of  finite  difference,  method  of  integral  relation  and  method  of 
lines.  Applied  to  smooth  bodies,  all  these  methods  can  have  satisfactory- 
results.  Only  because  the  nonstationary  method  must  take  steady  process  for 
time,  it  has  to  consume  a  great  deal  of  machine  time.  As  for  the  method  of 
finite  difference,  in  order  to  have  very  precise  result,  it  needs  quite  a 
number  of  net  points  and  large  machine  storage  capacity,  and  it  uses  more 
computing  time.  Compared  with  these  conditions,  the  method  of  lines  has 
more  points  of  excellence.  For  instance,  its  computing  method  is  simple,. the 
storage  capacity  it  needs  is  small  and,  using  only  a  few  rays,  it  can  bring 
about  satisfactory  result.  This  article  is  intended  to  report  our  work  of 
using  the  method  of  lines  to  compute  supersonic  encircling  flow  of  blunt- 
nosed  body. 

We  use  the  method  of  lines  to  make  broad  computation  of  supersonic 


encircling  flow  of  blunt  nosed  body.  The  objects  we  computed  include 

ellipsoid  of  various  axial  ratio  and  disk-analogie  bodies.  The  range  of 
incoming  flow  M  number  is  1.5  to  infinity.  Under  the  condition  of  axial 
symmetry,  besides  the  frozen  gas  of  y  =  1.4,  we  have  also  computed  balanced 
and  unbalanced  air. 

In  addition,  we  ,too,  use  the  method  of  lines  to  compute  the  flow  in 
supersonic  zone  and  the  pointed  conical  encircling  flow  with  attack  anele. 

To  the  results  of  computation,  we  make  multi-way  check,  and  all  show 
that  the  results  of  computation  by  using  method  of  lines  are  considerably 
satisfactory. 


2.  Tha  Way  of  Relaying  Questions 
2.1  Fundamental  Equation 

To  consider  the  inviscid  and  non-’neat  conducting  air  flow.  The  equation 
of  aerodynamics  in  spherical  coordinate  system  (r,  0,  cp)  is: 
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Here 


+ _ !_JL 

r  dd  r  dn  6  dip 


,  *,  »,  w  are  component  when 


di  dr 

speed  is  alona .  t ,  0,  9 direction,  p  is  pressure,  p  is  density  and  c  is 
speed.  In  the  equation,  all  quanity  are  dimensionless  quantity  and  their 
dimension  factor  are  respectively 

P  ~  P-Vi,  P  ~~  pm,  U,  *>,  W  ~  Vrn,  r  ~  R, 


Here  00  indicates  the  quanity  of  incoming  flow  and  R,  is  the  curvature 
radius  at  the  top  of  the  subject. 


Pbr  the  convenience  of  computing,  we  introduce  the  following 
transformation  of  coordinate, 


{ 


r  —  G( 0,  <p) 

He,  <p)-  g(s,  *>)’ 


0,  <p 


<p 


Here  r  —  c(3,  <p)  —  F(3,  <p)  are  equations  respectively  of  the  object 

surface  and  shock  wave.  Obviously,  in  (£»  v)  coordinate  system,  the  shock 
wave  is  in  the  plane  of  f  —  l  and  the  surface  of  the  object  is  in  the  plane 
of  S  —  0  .  Let 


-  (C.  +  $e»),  ? - (C,  +  e  -  F  -  G 

SID  0 


Due  to 


a  w  1  a  a 

Sr  e  8f  ’  dd 

1  a  -  f  a  , 

sin  3  dip  £  Q  |  sin  3  dip 

so  equation  (2.1)  can  be  rewritten  into 
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Or  written  into  solved  form 
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In  the  equation  , 
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2.2  Boundary  Condition 
(1)  Condition  of  shook  wave  On  the 
shock  wave,  the  shock  wave  relation  equation 
must  be  satisfied. 


P  +  pl'i  -  p»  4-  pmVi„ 

* - (l  - 

’  -  ?)"*•'- 
-  [i  -  eAniv„m 


(2.4) 
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Here  the  quantity  indicated  by  an  infinity  mark  is  wave-front  quantity  and 

that  not  indicated  by  an  infinity  mark  is  wave-back  quantity,  V.  is  the 
speed 

projection  o!^  along  the  direction  of  wave  normal  line,  h  is  han, 

are  direction  cosin  respectively  along  the  shock  wave  normal  line,  namely 


(2)  Condition  of  the  object  surface.'  On  the  surface  of  the 

object,  it  must  satisfy  the  condtion  that  the  normal  direction  speed  is  zero, 
namely 

q  uC  —  vC$  —  w  C*  —  0  (2.5) 
sin  6 

3.  'luraerical  Solution 

In  order  to  make  numerical  solution,  we  introduce  some  rays  to  the 
solution  zone.  For  instance,  at  0  direction  we  introduce  coordinate  surface 
of  0  =*  6,  —  ynptt  and  at  <P  direction,  we  introduce  coordinate  surface  of 

<p  mm  const  ,  then  we  take  the  intersecting  lines  of  these  coordinate 
surfaces  as  rays.  J bri>:  we  use  the  value  of  blow  parametre  on  the  ray  as 
nodal  point  value  to  construct  interpolating  polynomial  equation  and  then  to 


determine  the  partial  derivative  of  6  &  V  correspondingly.  Let  equation 
(2.3)  be  formulated  y.i  the  ray,  then  we  have  a  constant  differential 
equation  group  and  the  problem  of  marginal  value  will  correspondingly 
become  marginal-value  problem  of  constant  differential  equation  group.  For 
this  reason,  we  change  marginal-value  problem  into  initial-value  problem, 
and  work  out  solution  through  iteration,  namely  we  first  assume  that  the 
form  of  shock  ■wave  has  been  known, and, then  from  shock  wave  condition  (2.4), 
we  have  the  flow  parametre  of  wave-back.  Taking  this  as  initial  value  of 
integral  constant  differential  equation  group,  then  we  check  whether  the 
flow  parametre  of  the  object  surface  can  satisfy  the  condition  of  object 
surface  (2.5).  If  not,  we  adjust  the  shape  of  shock  wave  till  the  condition 
of  object  surface  is  satisfied.  The  integration  of  constant  differential 
equation  can  use  general  method,  such  as  quar trie-valence  Runge-Kutta 
method.  In  the  following,  we  shall  describe  some  specific  treatment. 


3,1  Equation  on  Axis  0  —  0 

Assuming  that  flow  field  is  symmetrical  with  V  “  °»*  plane  and  that 
axis  0  •» 0  is  always  in  a  symmetrical  plane.  From  equation  (2.2),  it  can  be 
seen  that  due  to  the  fact  that  ri a  6  appears  in  the  denominator,  the  equation 
at  6  -  0  must  be  given  a  treatment.  Let  *•-*<«»  °»  0)»  clearly  »(«,  0.  <P>  “ 
**(f)co*9>  “'({>  °»  ,,) --**(!)  iin  <p,  «CS,  °>»  0,0), 

p(f»  0.  ?)—  °»  0).  By  applying  Low-bi-ta  (transliteration  of  Chinese 

sound  and  it  may  be  a  Chinese  transliteration  of  Robert)  method  to  o/o  which 
appears  in  equation  (2.2),  we  can  have  the  equation  on  0  —  0  .  in  principle, 
it  will  do  by  taking  any  equation  from  <P  surface  randomly.  But  because  the 


6 


computing  error  of  numerical  value,  of  different  <p,  there  will  be  different 
results.  In  order  to  eliminate  such  incongruity,  we  make  integration  of 
those  equations  of  <P  from  o  -  *  to  induce  the  necessary  equation.  For 
instance,  we  use  cos q>  to  multiply  the  third  equation  of  equation  (2.2), 
and  take  off  the  fourth  equation  and  '  use  <P  to  multiply  it,  then  we  make 
integration.  Due  to 

(  a  (^-  cosip  —  sin  y'Wy  —  (  (ru  +  av  +  ^<p  t*  -~~ 

)•  \d£  dg  J  J»  dg.  dg 

—  (  (ocotro  —  «,sin  <p) 

Pi*  dg  P  dg 

(F,coi<j>  —  F.sinqp)^  —  —  ~{j,  [*'*  (  COI>’>  —  ^  <P coi  v) 
4-  —  ^2-  co»<pl  J<p  4-  —  in’*}  ■“  F* 

P  dd  J  2  i 


So  we  can  have 


And  analogously  we  can  have 


„*  Pal  +  «!  $£.  -  f t 

dg  P  dg 


PlL  +  JL$L  -  F? 
dg  p  dg 


„•$£.  +  p(PlL  +  a»  -  FT 

dg  \dg  dg  ) 

ap\ _  F. 

\dg  dgJ 

To  write  into  solved  form  of  |f,  •*•,  —,  we  can  have  the  necessary  equation 

Os 


p£_  _  CJ[  Fra*  -  p(rF?  4-  g*F*)1  4-  a*FT ' 
dg  a*3  -  r*V 

0fi.-.i[_£L+  0£_1 

dg  (}  I  a*  dg  J 

du  _  1  [  _*  r  dp  1  I 


9£.  =  _L[_£L+  0l1 
dg  c>  l  a*  dg  J 

du  1  f  r*  r  Q/>  1 

e7  7*  lFj  ”  ~Pdg\ 

du*  _  sh 

«*  1  P  dg  J 


In  it, 


a*-  i.  fa 


cos  ipdcp 


a  “=  ru  +•  aV* 
r*J  -  r'  -t-  a*1 


F*  - 
F?  - 
F*  - 
F*  - 


“  ¥  {"*  +  pJ."i-^  +  o^} 

~^{v*  LSr~**-rw} 

”  ^  ft  t"*  (lr  cosV  ”  I? s!n  ’>cos*)  +  7  |J-  «*'/>]'* 


3.2  Interpolating  Multinomial  Equation 

Because  the  flow  field  is  assumed  to  be  symmetrical  with  <P  “  °»  *  plane, 
it  is  only  necessary  to  have  solution  between  0  <  <p  <  * .  Between  0  *,  we 
introduce  4  +  1  planes.  And  at  the  sa~e  tine,  we  make  n  +  1  conical  surfaces, 
0  -  ek  -  comt  id,  -  0),  {  -  0,  1,  . . . ,  n.  They  intersect  wi  th  half  plane  <p  -  <p, 
and  ?>  —  *>.  +  *  to  form  (2n  +  1)  rays.  Noticin'*  that  the  flow  parametre  on 
the  plow  symmetry,  <p  —  «■— <p,  and  the  flow  parametre  on  <?  —  ?>,  +  *  are  equal  or 
different  by  one  symbol,  for  the  fixed  S,  we  can  utilize  the  value  of 
2n  +  1  rays  on  <p  ■“  <p,  &  <p  —  *  —  sPi  to  construct  2nth  order  interpolating 
multinomial  equation  Of  0. 

i  -  S  (2.7) 

i  «0 

g  indicates  flow  parametre.  For  the  purpose  of  savin'-  time  in  the  process 
of  computing,  we  do  not  first  compute  coefricient  oj.but  use  the  following 
computing  methods  instead.  Because  interpolating  function  and  its  derivate 

3 


■  "W  ■ 


can  be  expressed  as  a  linear  combination  of  function  values  on  interpolating 
nodal  point,  and  the  linear  combination  coefficient  is  only  related  with 
the  position  of  interpolating  nodal  point  and  the  position  of  interpolating 
point,  so  when  t’^e  nodal  point  and  interpolating  point  are  ^ iven ,  these 
coefficients  can  be  determined.  When  the  function  of  each  point  and  the 
derive he  value  are  computed,  it  will  do  to  use  these  coefficients  °r.d  the 
function  value  on  nodal  point  to  make  point  product.  Taking  equation 
(2.7)  as  example.  If  nodal  ooint  is  6m (m  —  0,  l,  •••,  2 »),  £(g)  at  0  of  some 
interpolating  point  will  be  computed.  Because  there  is  condition  at  nodal 
point, 

2  “  t-  (*»  “  0,  1 ,  •  •  • ,  2») 

i  »0 

gm  “  g(0m).  or  written  into  matrix  form: 

M  a  —  g 


Here 


M 


1 1  e,  oi 
1 1  e. . 


Therefore  we  have 

a  -  M-g 

Then  we  can  rewrite  the  equation  of  fW, 

g(0)  -  d*  ■  • 


o 


Into 

tie)  -  d*  •  (w-g> 

-  (m— d)*  ■  g 

—  b*  •  g  (2.8) 

Here  b  —  M~'*d,  gd*  —  (1 »  0>  *  *  * »  Evidently,  when  nodal  point  and 

interpolating  point  are  given,  M*  and  d  can  be  determined  and  then  we  can 
have  b.  Sinilarlv  because, 

S»W  —  ai*W~' 

»-0 

-  dr  •  • 

-  (M  — d,)*  ■  t 

Here  d,*  <->  (0,  1,  28,  •  •  • ,  represents  derivate  of  0,  so  to  compute 

derivate  can  be  of  an  analogous  treatment. 


So  far  as  q>  is  concerned,  when  0  .is  fixed,  we  can  utilize  the  value 
at  its  intersecting  line  with  Jfc  +  1  planes  to  construct  trigonometric 
interpolating  multinomial  equation  Of  ?*.  In  computing,  the  method  mentioned 

above  can  also  be  used.  Fbr  even  function,  we  take, 

> 

g  —  2 

• »  o 


Then 

TIow  there  is 


M*  - 


~  (w—d.)*  •  * 

i  i  —  i 
cot  <p„  cos •••  eosg>t  I 


<cos*f>g  «»*$>,•••  eo  s*<pj 
d|*  —  (0,  —  tin  qp,  •••,  — ^  tin  qocos  *_1qp) 
g*  “  (f«»  fi>  ”•» 


i 


10 


Fbr  -L  ytj<p,  —{  gconpJip  ,  3inil.Tr  exnression  can  also  be  written.lt  is 
not  necrssary  to  exanplify  them  here.  Fbr  odd  function,  we  take, 


Then 

Now 


*  - 1 

g  “  2  cos'?) sin  ?> 
l-Q 


(a/— d.r* 


A/* 


sin  <7>, 
sin  tp,  cos  <pi 


tin  ?>,• 


sin  q>, cos'  'ip,  sin  ?>,cos*_1?>,  •••  sin  ?>,-, cos' 


tin  ?>*- 


d,*  «*  (cosqp,  —  sin1?)  +  cor1  —  2)cos*"’?>*>nV  +  cos*-,?>) 


In  addition  to  the  methods  mentioned  above,  we  also  use  the  following 
methods  to  construct  interpolating  multinomial  equations.  Fbr  even  function 
of  <P  ,  we  take 

•  < 

g(£)  “  E  S  g.X‘)s‘cos',P 

•  =  0  1=0 

Fbr  odd  function,  we  take 

“'(I)  =  (2  S  WiXD&coi'fp)  simp 
'»  *0  >  =0  ' 

In  order  to  make  the  function  and  the  deriva+e  of  ®  at  0=0  in  some  sense 
be  single  value,  some  proper  condition  must  be  added  to  the  multinomial 
equation.  Now  we  try  to  describe  such  conditions. 

For  instance,  for  shock  wa^e  form,  • ,  _  f(fi,  we  naturally  require, 
when  0  —  o  j  it  has  no  relationship  with  This  means  that  we  require  when 
i 0  >  F*  —  0.  And  at  the  same  time,  in  order  to  warrant  that  the  normal 


line  has  definite  direction  and  when  8  ~  o  ,  the  normal  line  direction 
spherical  coordinate  should  have  such  a  form:  (fl>  bcosip,  —btio <p>  an<^  ^  ^ 
a  and  b  are  constants.  This  means  to  require  F„ -■  0  (j  **  O.  So  for  F,  the 
multinomial  equation  should  be: 

•  » 

F(0,  <?>)■“*  F «  +  F„0cos<p  +  2  2  Fii&cos'V  (2.9) 

.  «I  /  “0 

Obviously,  u,p,  and  p  should  also  have  the  sane  form  as  F.  The  above 
conclusion  and  the  form  of  v  and  w  can  be  obtained  as  well  in  the  following 
way.  Let 


g  “  foo(f)  +22  g‘! (S)&COt,9 

i  •  J  /  *  0 
W  k 

V  "*  »tn(f)cos tp  +  2  ‘>,Xi)Sicot’q> 

i  -  I  /  *0 

«’  ■=  —  <\,,(?)sin<p  +  (2  2  "'i/(Z)8icodq>)na  <p 

1  y-o  • 

Here  g  can  be  used  to  express  p,/>,u,F,  and  utilizes  the  property  o,  <p)  — 
4  cosy,  <«■(£,  o,  <p)  —  —  Wsin  y,  then  there  is  —  0  (/  =v  l),  wt,  «=  o(/  0),  vn  *=  —  «'*. 

To  combine  the  above  equation  with  equation  (2.2), we  notice  that, 

G(0,  y)  —  C,(0,  0)cosy, 

3-1  ■“  —  C#( 0,  0)siay 

sin  0 1 8—1 

(Here  the  object  body  is  symmetrical  with  <P  -  0,  *•  surface),  so  when  0 -*•<), 


we  have 


r°PK~zf  ~~  P»lG»o  +  5(Fu  "**  Fj 

0  in*  +  lsi£i£  -  \tm  F, 
dt  p»  dZ 


In  it 


f«0  ^  -  gj"  +  *(fn  -  dfiA 

Ae  Jt ) 


COS  tp 


-  (*e  ^Si  -  go;  ±  K £u 
V  '»c 


^iea\ 


sin  ip 


lim  F, 


*•  lim  Fa 

0~Q 


lit*.  _  f«4aA 
V  ^  di ) 


lim  r, 
•  -0 


«C  —  H»r,  *'0.  (  Cfc  +  f(f  „  —  G*,)) 

r»  ~  C(0,  0)  +  £IF„-  G(  0,  0)] 
C^~  ~GC<i,0) 


cl  — 


LP.& 

Poo 


To  write  the  ri~ht  end  of  the  above  equation  into  multinomial  equation  of 
cos  <p  (or,  in  addition,  to  multiply  it  by  sin  7)),  and  to  use  the  linear 
independence  of  1,  cos<p,  •  •  •  ,cos*t>  ,  ve  Can  reason  out  that  g  should  have  the 
form  of  equation  (2.9)  and  v  and  w  should  take, 

"  ♦ 

V  **  VH  COS  7>  +  V[0Q  +  l>)fi£o$‘<p  +  V]  2j  I'li^'cos’ip 

'  •  t  I  *  u 


«'  “  —  t/0,  sin  7)  —  t/,,6) cos  <i>  sin  7)  -+-  f  i<'i,0'cos'<p)  sin  q> 

'1  .1  i  =  o  ' 

At  the  same  time,  we  have  equation  on  Q  —  Oi 

‘0~?  +  r,pM~f  -  PvlCeo  +  f(F„  - 

“  —  (Poo  —  C(0,  0))lt'o1Pn  +  Pou(«'i:  +  2v„  +  2 «M)] 

P 

«o  —  “  ~  (Poo  —  C(0,0))t'o,(»u  —  f\>i) 

</?  Poo  *1 

a  rft,ot  _  Can  -f-  f  (  PM  —  Cm)  t/poo 

</!  P*  rf? 

—  —  (Foo  —  C( 0,  0))  [n.Ct'io  4-  r„  4-  #»)  + 

1  PiOJ 

•. “  rllf)  “  “  (F»  -  c(0* °»  •'«  (ft*  -  ^  *.) 


(2.10> 
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3.3  Iteration  Method 

As  what  has  been  stated  at  the  be^innin"  of  this  section,  we  shall  try 


to  solve  the  problem  of  mar~innl  value  through  iteration.  In  fact,  it  can 
be  interpreted  as  a  problem  of  solving  a  transcendental  equation  --roup .  That 
is  to  select  one  vroup  of  5b  (to  indicate  the  shock  wave  form  on  ith  ray). 

To  make, 

q,(F ,,  (;  “  1,  2,  ra) 

Here  q,  «*  o  is  the  boundary  condition  (2.5).  We  use  Newton  method,  namely  to 
use  alteration  quantity  8F ,  of  Fi  to  satisfy  the  following  equation: 

S'  J^SF,  -  (i-1. 

OFi 

partial  a 

Usually  there  is  no  way  to  express ^derivate  by  using  analytic  equation, 

so  we  use  numerical  value  method,  namely 

pq,  „  q,(.Ft,  •  •  -  ,  f,-i,  F,  -f-  AF,,  F..,,,  •  •  •,  F„)  —  S,(Ft,  •  •  •,  F„) 
dF,  AF. 

3y  using  this  method,  it  needs  m  +  1  times  of  integration  for  each  iteration, 

so  it  consumes  a  'Teat  deal  of  machine  time.  For  the  purpose  of  saving  time, 

we  can  use  the  simplified  Newton  method,  but  because  of  the  lack  of  accuracy 

in  msot  of  ,  the  speed  of  convergence,  therefore, can  possibly  become 

slow  In  order  to  make  r^L  more  accurate  without  increasing  much  of  the 

dF, 

volume  of  computation,  we  su~gest  a  method  as  follows.  Let  Q  indicate  the 
vector  constructed  from  q„  ••*»*-  R  is  the  vector  constructed  from  Fq, 

...,  5h.  Q  *  ^derivative  index  of  Q  to  R.  If  Q(Ro)  =  Qo»  ot^er 

of  Ri,  ...,  Ba,  close  to  Ro,  Q(%)  =  Qi(i  =  1,  ..-,»)  ***  been  known,  and 
if  _  R^(i  =  1,  2,  . ..,  m)  is  linear  independent,  then  Q’0  can  be  decided 


approximately  by  %  and  <&.  In  fact,  because 


So  there  is 
and  then  there  is 


P.  «=P.  +  Q',(R.  -  Rj  («  -  I,  2,  «) 

(Pi  Po,  *  *  *  »  Qm  P.)  Ob(  /?,  /?oi  •  •  • ,  Rm  —  Rf) 


Pi**  (P.  -  Po.  •••,  Qm  -  P.X*,  - 

thus  we  can  have  an  iteration  formula 

“  R.-(.R.-m  -  R„  -  R.xp._-  p.,  •  •  • ,  P.-,  -  P.)-'P. 

» —  in  +  i,  •••  (2 m 


Evidently,  using  the  above  equation  to  make  iteration  and  begin  with  Rj, 
%+2  to  solve  Pi.  •••»  Qm+i  needs  m  +  1  times  of  integration.  But 
thereafter,  one  iteration  needs  only  one  time  of  integration. 


Now  we  try  to  make  a  simple  estimation  of  the  speed  of  convergence. 
Obviously,  when 

Rm+l  s  (  Rl  Rm+l,  ■  *  ■  .  Rm  Rf»+l)  •“  Af  ■  E 

equation  (2.11)  becomes  difference  half  Newton  formula.  Here  E  is  unit 
matrix  and  is  pure  quanty.  Because  now  there  is, 

Qm+l  =  (Pi  —  P-+1.  Qm  —  P-+l)  —  P-+lR-+l  +  *) 

Then  there  is 

0™+',  -  R.,*,P-V,  +  o(AF!)P-*,  -  R«+iP»+i  +  o(^F) 

And  thereupon  we  can  have  an  estimation  equation  for  difference  half  Newton 
formula. 

R  Rm*l  —  R  R«i+I  t  ^m+lP'B  +  lPm+l 

“  R*  R«i+i  +  P*+tP»i+i  t  (^<»+iP«.+i  —  P»*i)Pi,+i 
-  o(  ||R*  -  R„+1||')  +  o(AF[|p-+.||) 

IIR*  -  R.ll  ^  ^IIR*  ~  R_ill‘  +  BAF|ip._,||  (»-  1,  2,  •••) 


or  written  into 


Here  R*  is  true  solution,  ||  II  indicates  mode,  and  A  and  B  are  suitable 
constants.  Similarly,  for  difference  simplified  Newton  method,  there  is  an 
estimation  equation, 

liR*  -  RJ  <  A[\ R*  -  R._,|)  +  £AF||£>._,||  +  C|!R.+1  -  (Ip.-, || 

(n  m  +  2,  •••) 

For  the  method  su-gested  by  us  there  is  an  estimation  equation, 

|| r*  — .R„||  <  <4|| R*  -  R.-t!l  +  D!iG»;iill  ,  *"P  "* 

Here  A,  B,  C,  and  D  are  constants. 

It  is  easy  to  see  that  when  A F  &  ||R*  —  R !|  are  of  same  quantity 
level  or  smaller,  difference  Newton  method  basically  maintains  the  speed  of 
convergence  of  Newton  method.  But  ||R»+,  —  R.-ill  is  generally  increased  as  n 
is  increased,  so  of  simplified  Newton  method  the  speed  of  convergence  is 
slow.  Because  fP™-ill  *«p  l|R.-i  ~  R.-t-J  —  «0),and  »UP  llR.-i  -  R.-i-.!l  is 

IO<«  !<<<• 

reduced  as  n  is  increased,  so  the  method  suggested  by  us  can  possibly 
converge  faster  than  simplified  Newton  method.  This  has  been  proved  in 
practical  computation. 

3.4  Selection  of  Initial  Shock  Wave  and  Interpolation  of  Object 
Surface  Quantity 

When  the  methods  mention&d  above  are  used  to  make  iteration,  the  success 
in  computing  will  depend  on  how  well  the  selection  of  initial  shock  wave  is 
made  .  Fbr  this  reason,  we  use  the  ready  results  according  to  the  way  of 
some  parametre  gradual  transition.  For  instance,  when  we  want  to  compute 


Mil 


the  flow  result  under  certain  attack  angle  •£*,  we  select  attack  angle  T)  as 
parametre.  After  the  result  of  has  been  obtained  (V  example  ^  =  0°), 
the  result  can  be  used  as  initial  value  to  compute  •  ’<• +  .  After  we  have 

had  the  result  of  rtU  +  Ajj  f  we  use  the  results  of  q,,  17,  4-  A17  to  obtain  initial 
value  of  j?o  +  ^4)  i)  br  way  of  linear  interpolation.  In  the  same  fashion,  till 
we  have  the  result  of  p*.  As  for  the  initial  shock  wave  form  which  is  needed 
in  the  berdnninr  of  computing,  it  can  be  secured  by  utilizing  the  result 
available  currently. 


Prom  equation  (2.3),  it  can  be  understood  that  on  the  object  surface 

a  =  0,  so  integration  cannot  reach  the  object  surface.  In  order  to  have 

the  quantity  of  object  surface,  we  use  extrapolation  method.  When 

integration  reaches  a  Certain  £*  (for  example  {•  —  0.1),  then  we  use 

the  £  values  of  a  few  neighbouring  points  to  extrapolate  the 
object  surface  quantity  for  example  using  the  values  of 
—  0.3, 0.2, 0.1  makes  a  quadratic  interpolation.  Here  we  would  like 
to  make  a  random  suggestion  that  if  the  other  computing  form, 
such  as  implicit  form  integration,  extrapolation  can  be 
completely  avoided.  For  a  situation  of  axial  symmetry,  we 
designed  another  form  to  compute,  and  the  result  proves  that 
it  is  a  success.  Here  we  have  no  plan  to  give  its  details. 


4.  Computation  Results 

We  have  made  broad  computation  on  blunt-nosed  body,  axial  symmetry  and 
three-dimension  space  encirling  flows  by  using  the  methods  mentioned  above. 
The  object  forms  we  computed  include  ellipsoid  of  various  axial  ratio  and 
objects  analogous  to  disk  (object  expressed  by  equation  *•  +  (**  +  y')*'1  -  i', 


»  >  2).  The  ran^e  of  incoming  flow  M  is  1.5  <A/„<ao.  For  the  situation 


of  axial  symmetry,  besides  the  frozen  "as  of  y  =  1.4,  we  compute  the 
balanced  air  as  well  as  the  unbalanced.  The  method  we  used  is  borrowed 
from  article  12  in  t’-'e  bibliography  appended  to  this  article.  The  patterns 
of  unbalanced  air  are  ^resented  in  another  article  of  ours.  The  precision 
of  our  computation  results  have  been  checked  by  several  liferent  ways.  One 
of  the  checks  is  made  in  computing  as  it  is  in  process.  For  instance,  we 
use  different  number  of  rays  and  different  integral  step  length  to  check 


the  relations  which  should  be  satisfied 
constancy.  Another  way  is  to  compare 
With  ot’-'er  results  acquired  from 
experiments  and  other  methods, such 
as  internal  relation  method.  All  the 
checks  we  made  indicate  that  the 
precision  of  our  computation  results 
is  satisfactory.  In  the  following  we 
shall  present  a  part  of  our  computation 
results . 

Figure  1  illustrates  the  forms 
of  shock  wave  and  sonic  line  of 
frozen  "as  spherical  flow  under 


by  flow  field,  such  as  maintaining 


^i'-ure  1  Forms  of  shock  wave 
and  sonic  line  of  frozen  "as 
encirling  ''low 


different  A/-  number.  From  the  Figure,  it  can  be  understood  that  the  forms 
of  sonic  line  are  of  two  different  types.  When  M  >  3,  the  limit  characteri¬ 
stic  line  is  the  second  family  characteristic  line  that  can  reach  sonic 


point  of  the  object  surface.  But  when  K^3,  the  limit  characteristic  line 
is  composed  of  the  first  family  characteristic  line  (which  enmes  from 
object  surface)  in  contact  with  sonic  line  and  the  second  family  character¬ 
istic  line  (  which  comes  from  shock  wave). 


Figure  2  Distribution  of  pressure 
along  object  surface 

Figure  2  shows  the  distribution 
of  object  surface  pressure  of  frozen 
gas  spherical  flow. 


Figure  3  Shock  wave  position 
and  sonic  line  form 

(M-  -  3) 


Figure  3  illustrates  Mm~  3,  r  “  1-4,  the  shock  wave  form  and  sonic 
line  position  of  different  objects.  For  very  blunt  body,  if  n  >2C,  shock 
wave  position  will  basically  maintain  unchanged.  After  the  contraction  of 
the  curvature  radius  of  object  surface  adjacent  to  sonic  point,  for  solid 
Mm  number,  beginning  with  a  certain  curvature,  there  will  be  torsional 
point  on  the  sonic  line. 


Figure  4  shows  shock  wave  position  and  sonic  line  form  of  balanced  air 
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spherical  flow.  The  conditions  of  incoming  flow  are  M  .  4 ,  p.  “  0  87  x  10> 

dyne/cm2,  p-  -  0.95  x  io_‘  g  /  en  .*&  M„  -  20,  p-  -  0.122  x  10‘  dyne/cm2, 
p»  --  0.192  x  io-J  g/cm3.  From  the  Figure,  it  can  be  seen  that  ionization 


Figure  4  Form  of  shock  wave 
and  sonic  line  of 
balanced  air  snherical 
flow 


Figure  5  Form  of  shock  wave 
and  sonic  line  of 
imbalanced  air 
spherical  flow 
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Figure  6  The  detatchment  distance  of  stationary  point  shock  wave 
following  m-  to  charge 


makes  the  situation  of  shock  wave  position  and  sonic  line  with  frozen  gas 

Y  —  1  4 

change  remarkably.  Shock  wave  moves  much  closer  to  the  object 
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surface. 


Figure  5  shows  shock  wave  position  and  sonic  line  forn  of  unbalanced 
air  spherical  flow.  The  conditions  of  incoming  **low  are  />-  —  0.947 x  10* 
dyne/cm^,  pm  —  0.123 x  10~5  r^/cm3t  R  =  5  cm.  '//hat  is  worth  of  attention  is  the 
special  forn  of  M«“  20  s^nic  line  at  the  place  of  shock  wave. 

Figure  6  shows  the  relationship  between  the  detatchnent  distance  of 
stationary  point  and  !il-  number.  For  frozen  flow,  when  M >  10 ,  it  remains 
unchanged.  For  balanced  air,  following  the  increase  of  Mm  ,  the  change 
of  detabchment  distance  apperas  to  be  not  unique. 

Using  5  rays  to  compute  ancirling  floy  of  axial  symmetry. 

Figure  7  shows  that  of  the  ellipsoidal  flow  of  s  “  *-5  when  *v/-  —  3  &  4  , 
shock  wave  and  sonic  line  in  symmetrical  plane  will  follow  the  change  of 
attack  angle  *}.  Figure  9  shows  that  the  object  surface  pressure  in  symmetric 
plane  will  follow  the  change  of  attack  angle.  Also  Figure  n  shows  position 
of  stationary  point  and, in  accuracv,  stationary  angles  of  Mm  —  3  &  A/»  —  A 
are  overlapping.  When  attack  angle  is  chanr-inm,  it  moves  alon^  object 
surface  by  almost  the  same  speed.  Figure  S  shows  that  shock  wave  and 
sonic  line  in  <*>■”*/ 2  plane  follow  the  change  of  attack  angle.  In  the  Figure, 
it  can  be  seen  that  the  chan~e  of  shock  wave  form  is  slow  and  the  change  of 
sonic  line  is  faster. 


Figure  10  shows  the  encirlin^  flow  of  disk-analogue  object  of 


Figure  7  Shock  wave  form  and 
sonic  line  in  sym¬ 
metrical  plane 
following  the  change 
of  attack  angle 
(S  -  1.5) 

(1.  sonic  line,  2,  shock  wave, 

3.  direction  of  incoming  flow, 

4.  stationary  point,  5. result 
from  article  12) 

from  each  <p  surface  we  take  4  rays, 
take  13  rays  altogether.  The  z  axis 
symmetrical  plane  of  flow  field,  and 


Figure  8  Shock  wave  form  and 
sonic  line  following 
the  change  of  attack 
angle 

(M_-  4,  S  -  1.5) 

Mm  “  5.8,  n  20  shock  wave  form 

and  sonic  line  position  on  symmetric 
plane  under  different  attack  angle. 
Figure  11  shows  the  distribution  of 
object  surface  pressure  in  symmetric 
plane . 

Wien  we  compute  space  encirling 
flow,  we  take  four  <P  surfaces  .and 
Because  the  axis  line  is  common,  we 
of  coordinate  system  is  placed  in  the 
,  ror  the  convenience  o*'  computing,  we 
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Figure  10  Shock  wave  form 
and  sonic  line 
position  on  sym¬ 
metrical  surface 

<*/.  -  5.8,  «-20> 


Figure  11  The  distribution  of  object 
surface  pressure  in 
symmetrical  surface 

(A/.- 5.8,  »-20) 


About  the  problem  of  whether  the  greatest  entropy  value  on  object 


r  * 


surface  can  be  reached,  from  our  experience  of  computing,  we  think  that  for 
a  certain  object  in  a  certain  error  range,  the  greatest  entropy  can  be 
reached. 

In  order  to  check  our  computation  results,  we  use  several  different 
kinds  of  ways.  For  instance,  for  axial  symmetrical  enc’rlin"  ''low,  we  use 
5  rays  and  3  rays  respectively  to  compute  and  the  result  shows  that  error  is 
no  more  than  1%.  We  also  use  different  integral  step  length,  for  example, 
for  spherical  flow  of  A/«  —  6,  r  —  1  < ,  from  shock  wave  to  object  surface  we 
integrate  10-step,  20-step  and  30-step.  The  relative  error  of  10-step  and 
20-step  is  no  more  than  0.35%  and  between  the  results  of  20-step  and  80- 
step  there  are  at  least  three  sane  effective  difrits.  This  means  that  we  do 
not  have  to  worry  about  the  increase  of  roundingoff  error. 

We  compare  the  results  of  — 3,  8  —  1.5,  i?  ■»  15°  with  those  of  Telenin^^ 
they  are  completely  identical  as  showed  in  Figure  1.  For  ellipsoidal  and 
spherical  frozrn  flow  of  '-5  ,  our  computation  results  have  three  coincide 
effective  dibits  with  the  results  Belotserkovskiy  obtained  bv  uslne  integral 
relation  method. 


We  also  made  integration  check  and  examine  the  accuracy  of  interval 


j  j  pp~'fju  •  do  •  0 

Here  x±  indicates  three  axial  unit  vectors  under  rectangular  coordinate, u 
is  the  projection  of  velocity  vector  u  at  three  directions,  a  is  a  curved 
surface  containing  no  shock  wave.  Under  the  condition  of  M*»  ••  4,  B  ■“  1-5, 
and  >?  “■  10°,  20e,  the  intonation  result  can  be  found  in  Table  1.  Besides 
the  total  intenation,  Table  1  also  b’-'ows  the  intenation  on  shock  wave, 


Table  1  Integration  Check  « .5, 4  «■  0.6$ 


*  Mass 

* 


Moss 

*  *%&&& 

Entropy 


Shock  wave 
surface 

. 

Object 

surface 

Conical 

surface 

Total 

-0. 64327x2 

0.00066x2 

0.64517x2 

0.00256X2 

0.04595X7 

-0.00731X2 

-0.03849x2 

0.00015  X2 

0 

0 

0 

0 

0.67039  X2 

-0.32282  x2 

-0.34862X2 

-0.0U105X2 

-0.05789x2 

0.00006X2 

0.05809x2 

0.00026x2 

-0.63088x2 

0.00084x2 

0.63298X2 

0.00294  x2 

0.09011X2 

-0.01579x2 

-0.07398x2 

0.00035X2 

0 

0 

0 

0 

0.65271X2 

-0.30987x2 

-0.34423x2 

-0.00139  X2 

-0.056379x2 

0.00008  x2 

0.05663X2 

0.00033  X2 

(■>  x,y,z  direction  momentum) 


object  surface  and  conical  surface.  PYom  Table  1,  it  can  be  seen  that  the 
computation  results  are  accurate,  and  there  is  error  only  at  the  third  digit 
of  the  integration  on  shock  wave  and  conical  surface.  The  object  surface 
condition  is  well  satisfied  and  it  is  at  lCP^  numerical  level. 

In  addition,  we  also  compute  the  total  energy  on  all  nodal  points  and 
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and  entropy  of  each  nodal  point  on  object  surface.  Prom  Table  2,  it  can  be 


Table  2  Error  of  Total  Energy 
The  accuracy  value  of  _  <  w  +  +  w> 

r  t  '  *  2 

+  is  0.65625 


V 

Total  energy 
maximum 
deviation 

’  Relative 
error 

0* 

0.0010 

0.15% 

5* 

0.0019 

0.30% 

10* 

0.0034 

0.51% 

15* 

0.0051 

0.78% 

20* 

0.0087 

1.33% 

seen  that  for  -  4,  6  —  1.5,  when 
n<  15°,  the  relative  error  of  total 
energy  is  less  than  1%.  Prom  Table  3> 
it  can  be  seen  that  for  A/„  —  4,  «  —  1.5, 
and  *?5£15°,  the  entropy  of  object 


Table  3  Ehtropy  Value  of  Nodal 
Point  on  Object  Sur¬ 
face  (Inp—  rlnp) 

When  «  -  i.5,  3  - 

o. 67.9-0,  the  accuracy 


0* 

10* 

20* 

0* 

0 

-2.31903 

-2.31537 

-2.30611 

11* 

0 

s 

-2.31853 

-2.31396 

-2.31974 

-2.31833 

-2.32091 

-2.30865 

-2.30962 

-2.31361 

-2.31858 

22*  i 

0 

i* 

* 

-2.31951 

-2.32010 

-2.31992 

-2.31939 

-2.32077 

-2.32559 

-2.32117 

-2.31759 

-2.32111 

33* 

0 

t- 

I* 

n 

-2.31723 

-2.3)276 

-2.31593 

-2.31851 

-2.32095 

-2.30283 

-2.31361 

-2.31903 

-2.32921 

surface  is  different  only  by  1  at  the  third  effective  digit,  and  for'*  20°» 


there  is  only  a  difference  by  3  at  the  third  digit. 


In  summary,  using  method  of  lines  to  compute  encirling  flow  of  smooth 
bodies  can  produce  very  satisfactory  results. 

Comrades  Rui  Wei-ming  participated  in  part  of  this  work.  He  Jiao-min 
gave  us  significant  help,  and  Feng  Kang  once  enthusiastically  led  us  to 
work  on  this  project.  Here  we  thank  them  all. 
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